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We develop a theoretical description of non-adiabatic Josephson dynamics in superconducting 
junctions containing low energy quasiparticles. Within this approach we investigate the effects 
of midgap states in junctions of unconventional d-wave superconductors. We identify a reentrance 
effect in the transition between thermal activation and macroscopic quantum tunneling, and connect 
this phenomenon to the experimental observations in Phys. Rev. Lett. 94, 087003 (2005). ft is 
also shown that nonlinear Josephson dynamics can be defined by resonant interaction with midgap 
states reminiscent to nonlinear optical phenomena in media of two-level atoms. 
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With the advent of superconducting qubits p~H2] a gen- 
eral interest has grown towards realization of macroscopic 
quantum dynamics in superconducting weak links. The 
superconducting qubits developed so far are based on 
Josephson tunnel junctions of conventional superconduc- 
tors. A conceptually interesting and practically impor- 
tant question is whether other types of Josephson weak 
links, such as junctions of high temperature supercon- 
ductors, and mesoscopic metallic or semiconducting weak 
links can be employed in qubit circuits. The central as- 
pect of this problem is to understand the role of low en- 
ergy electronic states usually present in such junctions. 
The low energy quasiparticles are driven away from equi- 
librium by temporal variation of the superconducting 
phase across the junction, and produce a non-adiabatic 
contribution to the Josephson current. This effect is com- 
monly considered to result in dissipation, and decoher- 
ence of qubit states. However, examples from nonlinear 
optics show that resonant interaction with localized elec- 
tronic states (two-level atoms) may generate a nonlinear 
dispersion of electromagnetic modes leading to a variety 
of nonlinear phenomena involving coherent energy ex- 
change between macroscopic and microscopic variables 
0]. This kind of nonlinear phenomena, whose origin dif- 
fers from the nonlinearity of the adiabatic Josephson po- 
tential, has never been studied in the context of macro- 
scopic Josephson dynamics. 

In this Letter we investigate the non-adiabatic Joseph- 
son dynamics in artificial grain boundary junctions of 
high temperature superconductors [5], which is caused 
by interaction with superconducting surface bound states 
(midgap states). The midgap states (MGS) situate at 
zero energy in the middle of the superconducting en- 
ergy gap [B], and are fundamentally connected to the 
unconventional d-wave superconducting order parameter 
in these materials pj [5] . We find that interaction with 
the MGS has implications in both the imaginary time dy- 
namics (tunneling) and the real time nonlinear dynamics 
of the junction. First, we show that the MGS are capa- 
ble of significantly affecting the transition between the 



thermal activation and macroscopic quantum tunneling 
(MQT) decay of Josephson current state inducing multi- 
ple, forward and backward, transitions between the two 
regimes. We suggest that such a reentrance phenomenon 
underlines the experimentally observed [5] anomaly of the 
switching current rates. Secondly, we show that the non- 
linear resonant response of d-wave junctions may be en- 
tirely caused by the nonlinear dynamics of the MGS, and 
lead to a bifurcation regime with an explosive growth of 
the response amplitude. These findings are made within 
the framework of a general theoretical description of the 
non-adiabatic Josephson dynamics in junctions contain- 
ing low energy quasiparticles, developed in this paper. 

The special role of the MGS is explained by their dis- 
crete energy spectrum, and pairwise coupling to the tem- 
poral variation of the superconducting phase. Tunnel- 
ing spectroscopy data [10] as well as observation of a 
7r-junction transition [TT] provide experimental evidence 
for the MGS existence. The equilibrium properties of 
MGS and their role in the dc Josephson effect are well 
studied in the literature (see reviews [H] H3] and refer- 
ences therein). The multiple degenerate zero energy level 
of the MGS splits into a narrow band under the effects 
of tunneling and anisotropy of the d-wave order parame- 
ter, A(kp) = Aocos(2(9). Due to the small bandwidth a 
thermal saturation of the MGS occurs at relatively low 
temperatures that may be comparable to the MQT tran- 
sition temperature. This saturation effect accompanied 
by the decrease of the MGS-induced dissipation under- 
lines, as we show, the reentrance effect in the MQT tran- 
sition. In junctions with atomically smooth interfaces, a 
large fraction of tunneling electron trajectories contains 
hybridized MGS pairs. The two-state Rabi dynamics and 
the MGS saturation at large driving amplitudes define 
the nonlinear property of real time Josephson dynamics. 

MQT transition temperature. We start with the discus- 
sion of the effect of MGS on the MQT transition temper- 
ature. We follow the method of Ref. [14], based on the 
analysis of the imaginary time dynamics of phase fluctu- 
ations, S(p(r), around the steady phase difference across 
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the junction, tp = at the top of the barrier of the 
tilted Josephson potential. In this method, the MQT 
transition is manifested by an instability of the phase 
fluctuations described with an effective euclidian action, 

SeffM « Seg[tpb] + J2n A(^, W n )^n^-n, V n = 2imT 

(ks = h = 1). The transition corresponds to the change 
of the sign of the kernel, A(tpb, iv\), and the temperature 
is given by the equation k(tpb,ivi) = 0. 

To derive the effective action for the superconducting 
phase, we consider the partition function of d-wave junc- 
tion, Z = J V(pV 2 ip e -(S v [ v ]+s.^[<pM) 7 an d perform in- 
tegration over fermionic variables ijj [13]. Here — 
J d,T[(C/8e 2 )ip 2 — I c ip/2e] is the macroscopic part of the 
action contributed by the charging energy of the junction 
capacitance, C, and the inductive energy of the biasing 
current, I e . Furthermore, = J dr J dr tp(d r + H + 
(i/4)sign(x)<£>)^ is the microscopic part of the action, as- 
sociated with the mean-field Hamiltonian of the super- 
conducting electrons, T-L, the last term provides electro- 
neutrality within the electrodes |16j . 

We perform the integration by choosing a general 
method suitable for all kinds of junctions regardless 
of their transparencies or presence of localized surface 
states. We separate the spatial problem from the tem- 
poral one by introducing a basis of instantaneous eigen- 
states of electronic Hamiltonian, Ufa — E-ifa, i/j(t,t) = 
'^2 i fa(r; ip)ai(r). The Fermionic action then becomes, 
5,/, = / drJ2ij o-iG^aj, where G" 1 = d T + Hij(tp, (p), is 
the inverse Green's function of the effective Hamiltonian, 



Hi. 



EiSij 



itpAij] (1) 
A i:j = [fa,id v fa) - (l/±)((/>i,sigD.(x)(T z <f>j) (2) 

is the matrix element of quasiparticle transitions induced 
by temporal variation of the phase. The effective action 
has the form, S e g[ip] — S v — SplnG -1 . 

The saddle point solution is given by equation, 5S e g = 
0. For the fermionic contribution we have, <5SplnG _1 = 
(l/2e)Sp (ijGStp), where 

i J = 2e(d <p E + i[E,A]] (3) 

is the Josephson current operator |16j . see Appendix. 
At the static saddle point, — G°(r, r) = h°(E) is the 
equilibrium density matrix commuting with E, there- 
fore only the diagonal (adiabatic) part of the current 
operator contributes to the Josephson current, Ij d {tp) = 
2eY^d v Eir^, that defines <p b , lf(<p b ) -I e = 0." 

The non-adiabatic effect is described by the 
second functional derivative of the fermionic ac- 
tion, (l/2e) 2 Sp (6<pijG°ijG 5<py and the fluc- 
tuation kernel reads (see Appendix), A(w n ) = 
(G/8e 2 ) (v 2 ~u} 2 - iUnjoliVn)). Here = 

— {2e/C)d v Ij d is the plasma frequency at the bar- 
rier, and 



,-|Ai| 2 (n?-n?) 



is the quasiparticle response; £ij = Ei —Ej, = np(Ei) 
is the Fermi filling factor, all functions are taken at if = 

Up to this point the derivation is general, and Eq. Q 
applies to all the quasiparticles. At small frequencies, 
however, only the MGS and itinerant nodal quasiparticles 
[17] are relevant. Furthermore, the MGS contribution has 
more pronounced temperature dependence compared to 
the nodal states because MGS have a small bandwidth, 
e m <C Ao- Focusing on the more interesting effect of 
the MGS, we truncate Eq. Q to the MGS subspace. 
The matrix elements, Aij, only connect MGS pairs of 
the same electronic trajectory while transitions among 
the trajectories are forbidden due to preserved transla- 
tional invariance. Parameterizing the MGS pairs with 
the angle, 9, between the incidental wave vector ]s.p of 
the respective trajectory and the interface normal (see 
top inset Fig 1), and denoting, e(6) — E\(6) — E 2 {9), 
A{0) = A12, we present the equation for the transition 
temperature on the form, 



V —U) h — 



\e 2 S 



eA 2 {n\ 



= 0, 



(5) 
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(4) 



where angle brackets indicate the average over the Fermi 
surface, S is the junction area, v = 2irT. 

The temperature dispersion of the MGS term in Eq. §5§ 
is primarily defined by the Fermi filling factors and the 
resonant denominator, while the particular form of the 
smooth functions e(6) and A(6) plays a secondary role. 
This allows us to formulate an analytical model equa- 
tion for the transition temperature, thus circumvent- 
ing the difficulty of evaluating anisotropy of the MGS, 
which generally can only be done numerically. By re- 
placing eA 2 {9){de/d6)^ 1 with a constant, we get Eq. ^ 
on the form, F(x) — e 2 n a; 2 (l + rjf(x)) — oj 2 = 0, where 
f(x) = Jq 1 dyta,iih('Ky/2x)(x 2 +y 2 )~ 1 , and a; = vje m \ 77 = 
8an J 'R n Ce m is the coupling strength, R n = n/e 2 S(D) is 
the normal junction resistance, and a ~ 1 is a geometry 
specific constant. The latter estimate is obtained from 
the scaling, e m oc \/DAo, and A oc y/~D, in the limit of 
small transparency, D< 1, extracted from the analytical 
equations for the MGS spectrum and transition matrix 
elements , see Appendix. The advantageous property of 
this analytical model is that it applies to junctions with 
interface faceting, which is taken into account by average 
values of the model parameters, 77, e m , and uib- 

Numerical solutions to the modeled Eq. ^ are pre- 
sented in the inset to Fig. [T] They demonstrate split- 
ting of a single critical point into three critical points 
manifesting the reentrance effect. The bifurcation of 
the solution to Eq. ^ occurs at the coupling strength, 
rj = 25, and the barrier frequency, oij = 3.45 e m . This 
phenomenon can be understood as a reentrance effect: 
At high temperature the thermal activation undergoes 
a transition to MQT in the absence of interaction with 



3 



MGS since the MGS are saturated; with lowering tem- 
perature, the MQT rate decreases because of increased 
interaction with MGS, and thermal activation takes over; 
then it undergoes the second transition to MQT in the 
presence of interaction. This finding constitutes the first 
main results of this paper. 

In the experiment with a tilt YBCO junction [9] an 
anomalous temperature dependence of the Josephson 
current decay rate has been observed, which can be inter- 
preted in terms of the reentrance effect: transition to the 
MQT regime at T x m 135 mK is interrupted, at T 2 « 90 
mK, by reentrance of the thermal activation, which then 
undergoes the second MQT transition at T 3 w 45 mK, 
as sketched on Fig. [T] To make a quantitative compari- 
son we fit the three experimental transition temperatures 
by adjusting the average model parameter values, r], e m , 
and u>b, see Appendix, as shown on Fig. [I] Including the 
stray LC-oscillator of the experimental setup [25] does 
not make any qualitative difference but rather insignif- 
icantly (within 20%) shifts the parameters values. The 
best fit is eventually achieved for the values, e m ps 320 
mK, ufy rj 1.7 K, ujp « 2.5 K, and C « 36fF, assuming 
experimental values of the critical current, 1q = 1.4/jA, 
and the switching current, I e 0.9/(7. Given the exper- 
imental junction transparency, D ~ 10~ 4 , we are able 
to evaluate the maximum energy gap at the interface, 
Ao ss 16 K. The geometrical constant in the equation for 
77 is estimated for the experimental value R n — 500 O, 
a « 1.5, as expected. 

In our discussion the temperature dependence of the 
adiabatic Josephson potential has been ignored. This de- 
pendence, also originating from the thermal saturation of 
the MGS band, may play a role in junctions with large 
capacitance where it may modify, as shown in [3], the 
thermally activated decay rate and provide an alterna- 
tive explanation to the experimentally observed feature. 

Consistency of our non- adiabatic reentrance scenario 
with the experimental observations strongly indicates in- 
volvement of the MGS pairs in the macroscopic dynamics 
of the junction. Moreover, it provides us with valuable 
information about the microscopic MGS parameters. 

Nonlinear resonance Josephson dynamics. To inves- 
tigate the real time Josephson dynamics, one needs to 
generalize our approach to non-equilibrium states. This 
is done by considering the partition function defined 
through the action on the real time Keldysh contour [22] . 
Then proceeding, as before, by introducing the instan- 
taneous basis, we derive the equation for the Keldysh- 
Green's functions, G ab , [id t - H(tp a , ip a )]G ab (t - t') = 
aS ah 5(t — t'), with the same Hamiltonian as in Eq. ([l]), 
here a, b = ± label the forward and backward branches 
of the Keldysh contour. The semiclassical dynamics of 
the superconducting phase is given by the least action 
principle, (6/5x)S g[(p, x] x =o = 0, formulated in terms of 
the Wigner variables, ip a = <p + ax/2 [53]. Calculating 
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FIG. 1. Reentrance effect in MQT. Sketch of temperature 
dependence of decay rate (wide shadow line) illustrates the 
effect featuring three transitions between thermal activation 
and MQT regimes. Experimental transition temperatures are 
given by zeros of function F(x), defined in the text (blue 
line) for n = 38. Lower inset shows development of non- 
monotonic feature of function F(x) with increasing n, at 77 > 
25. Upper inset illustrates junction geometry and scattering 
electron trajectory (dashed line). 

the functional derivative, we get, 

+ Tr (tjpj = I e , Ij = 2e(d lp E + i[E, A]). (6) 

Here p(t) = {^/2i)'}2 ia G aa (t 1 t) is the non-equilibrium 
single particle density matrix, which satisfies, by virtue 
of the equation for G ab , the Liouvillc equation, 



ip=[H,p], H = E-(pA. 



(7) 



Eqs. ([§ and Q are exact in the semiclassical limit, and 
give a general description of the non-adiabatic Joseph- 
son dynamics in all kinds of junctions. These equations 
constitute another main result of this paper. 

For the MGS pairs, Eq. |7|) reduces to the Bloch equa- 
tion for the two-level density matrix parameterized with 
the angle 6. In this case, Eqs. Q, §f§ become analogous 
to the ones describing electromagnetic modes in a cavity 
embedded in a medium of two- level atoms [4] . The most 
interesting is the case of the resonant excitation of the 
MGS pairs, which corresponds to the Josephson plasma 
frequency lying within the MGS band, lo p < e m . Suppose 
a small oscillating biasing current is applied to the junc- 
tion, I e cosujt, with frequency slightly detuned from the 
plasma frequency, 5 = uj — u) p -C oj. The resonant dynam- 
ics of the superconducting phase, <p{t) = Re^^e™ 4 "*), is 
described by the averaged equation for slow varying com- 
plex amplitude, 

- 2i(p u + [-25 + 7 (r)] (f u = eI e /uj p C, (8) 

where 7 = 7' + ij" is the nonlinear MGS response, 



7'M-7o + ^£^r 2 , 7 "(r): 



yJ(rAwY + r 2 



(9) 
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(the nonlinear adiabatic term is dropped from Eq. Q 
to emphasize the MGS effect). In Eq. Q the bar in- 
dicates the resonant values, r = \<p u \, and the quan- 
tity 7q refers to the linear MGS response given by the 
analytical continuation of Eq. Q to real frequencies, 
iv — > iv + iO. The response is calculated (see Appendix) 
by solving the Bloch equation ^ assuming the MGS adi- 
abatically following, in the rotating frame, the evolution 
of the phase amplitude, and adding small decoherence 
rates r 1; r 2 -C e m . The MGS decoherence is induced, 
e.g., by scattering to the itinerant nodal states by the 
facet edges or other rare inhomogeneities, leading to the 
MGS intrinsic broadening, Y — yTiT^- The dissipative 
part of the linear response is^estimated^ 

7o(w,T)~ ^^tanh— . (10) 
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It gives the frequency independent quality factor at zero 
temperature, Qmgs = w /7o ~ £ m RnC. It is instruc- 
tive to compare this result to the damping effect of the 
nodal quasiparticles, Q no dai ~ A R n C > Qmgs, see 
Appendix (cf. [MS])- 

Equation (|9| provides an extension of the linear re- 
sponse equation Q to the nonlinear region, when the 
Rabi frequency of MGS transitions exceeds the MGS in- 
trinsic width, rAuj > T. In this nonlinear regime relevant 
for narrow MGS levels the stationary response amplitude 
as function of detuning is defined by relation, 

S = ~ 7 '(r) ± ^{eIJCu p y - [ 7 "(r)rf. (11) 

The response demonstrates the bifurcation regime shown 
in Fig. [2] which is typical for nonlinear oscillators, 
but here is entirely controlled by MGS characteristics 
rather than adiabatic Josephson potential. The bifur- 
cation appears at very small driving currents, I e — 
(I e /2Ic)QMGS ~ T/ujA < 1. The most striking fea- 
ture of the response is the explosive growth of the peak 
amplitude, r max = 7 e [l - (Ie//*) 2 ]" 1 / 2 , for the driving 
current approaching the value I* = T/ujA. This effect 
is caused by the MGS saturation at large driving am- 
plitudes, which is manifested by decreasing damping in 
equation (|9| . The divergency is smeared by adding small 
damping, e.g., by nodal quasiparticles, and changes to 
a steep dependence asymptotically approaching the line, 
r ma x = (I e - I* e ){Qnod/QMGs)- The Rabi dynamics of 
the MGS should be more clearly exposed in the time re- 
solved experiments. 

In conclusion, we considered the effects of midgap 
states on Josephson dynamics in d-wave superconducting 
junctions. The analysis is based on the developed gen- 
eral theoretical framework for non-adiabatic Josephson 
dynamics in junctions containing low energy quasipar- 
ticles. We identified a reentrance effect in MQT, and 
connected that to the experimental observations. We 
also investigated the nonlinear dynamical response of the 
junction caused by coupling to nonlinear MGS dynam- 
ics. By analyzing the experiment [5] in terms of the 




FIG. 2. Effect of MGS on nonlinear resonance response of 
the junction. Response amplitude as function of detuning 
is shown for different amplitudes of driving current. Inset: 
maximum response amplitude as a function of driving current, 
dots indicate current values in the main figure. 

interaction with MGS, we found that the MGS band- 
width in the experimental junction is smaller than the 
Josephson plasma frequency, e m < uj p . This implies that 
the resonance condition for MGS excitation is not ful- 
filled, and MGS should not affect the real time junction 
dynamics, which thus would be similar to conventional 
Josephson oscillators. The quality factor is then defined 
by the nodal quasiparticles, and is estimated from the 
experimental data, Q no d ~ (Ao/w p ) 2 ~ 40. In order to 
increase this factor the strategy would be to increase the 
ratio, Ao/wp, which is, however, impossible beyond the 
limit, Qnod ~ provided MGS remain off-resonance 

(e m ~ a/DAo < u}p). Exceeding this limit necessarily 
implies resonant excitation of the MGS, and establishing 
the nonlinear regime described in this paper. 
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Appendix 

In this appendix we present details of the derivation of (a) effective action for superconducting phase, multiple 
critical temperatures for transitions between thermal activation and MQT regimes under interaction with MGS; (b) 
MGS energy spectrum and transition matrix elements, MGS linear response and comparison to the damping effect of 
nodal quasiparticles; and (c) the nonlinear junction dynamics under resonant interaction with MGS. 

Reentrance effect in MQT 

In this section we derive the dispersion equation for small phase fluctuation in imaginary time used to evaluate 
the crossover temperature between thermal activation and MQT decay of persistent Josephson current. To this end 
we shall also need to derive equations defining the MGS characteristics: energy dispersion equation, and interlevel 
matrix elements, and discuss MGS general properties, and present some explicit analytical equations. 

Imaginary time approach 

Starting from the imaginary time representation of partition function and performing integration over the fermionic 
fields, one obtains the equilibrium partition function as a path integral over the phase with an euclidian effective 
action S^P, 



SplnG" 1 , (12) 



C 



2 p 2 + U cxt (<p) 
8e z 



Z= J Vtpe-^M, Sf[<p} = J dr 
where f7 ex t(<^) = —(I e /2e)ip is an inductive energy of a biasing current I ei and 

6Z l = (d T + fly {if, tpj) , fly = 5ijEi - iipAij . (13) 



The matrix E in this equation is constructed with the eigen energies, £y — EiSij of the microscopic junction 
Hamiltonian, 

Ufa = E i( t> h (u) 

HV) 2 



u 



H + V 



a z + Ae lx a+ + Ae~ lx a- , (15) 



2m 

where A denotes a non-local operator, A0 = J dr' A(r, r')<p(r') = j dk A(k, r) j dH <fr(r + R)e lk R , and 

-f, rei? 

denotes the phase of the order parameter in the left (L) and right (R) electrodes; V(r) represents an interface potential. 
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The matrix of operator A is constructed with the matrix elements of the transition between the basis states, 
Eq. gg, 

Aij = (fa, [id v - (l/4)<r,sign(x)] fa). (17) 
An alternative representation of this matrix is given through the relation to the matrix Zy = (</>j, Ij (j>j), 

* j= 3, (18) 



iAi=-^ 
^ " 2ee tf : 



where £y = -Ej — Ej . As it is shown in the next section, this matrix represents the Josephson current flowing through 
the junction, and is connected to the current density operator via equations, 



Iij = / dn- jij(r), 
Js 



2m 



(zV-iVO^(r)^(r') 



(19) 



Current operator 



Here we present a proof for the interpretation of Ij as the quantum mechanical current operator. To do so we 
identify the current density operator as 



jij( r ) 



( 4 V- i V')0l(r)^(r') 



2m 



(20) 



The current through the interface, S, is given by 



Iij = I dn ■ jij(r). 
Js 



If we consider the system as an infinitely large loop we can use the fact that no current is flowing through any other 
part of the surface of the superconductor so we may extend the surface S around the whole superconductor and use 
Gauss law: 



21, 



f dyV.jiiW- f dVV-jiji*) 

JrEL JrER 



(21) 



From the explicit form of the BdG Hamiltonian ( 15 ) one finds the relations, 

-V 2 ^(r)]t^(r)-4(r)[-V 2 ^(r)] 



-*V) •j ii (r) 



e 
2m 



(Ei - E,)4(r)a z ^(r) + 4(r)[H, 



(22) 



The last term (the commutator between the " charge operator" a z , and the quasiparticlc Hamiltonian) can be rewritten 
as 



[a z ,H] = 2id x H = 

The current operator then becomes 

Iij = H)2e 

By differentiating the eigenvalue equation T-L^i = Ei<pi wrt ip one obtains the following identities: 



-4id v n, reR 



1 



,id v H<j)j) - (Ej - E l )-((j> i ,sign(x)a z <j) j ) 



i ,id v ,%4>i) = id^Ei, 

l ,id v M4> ] ) = (Ej - Ei)(<f>i,id 9 <t>j), i ^ j 



(23) 



(24) 



(25) 
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From this one sees that the current matrix elements are given by 

In = 2ed v Ei, 

or in explicit matrix form 

I = 2e(d v E + i[E,A\Y (27) 

Quasiclassical wave functions 

In this subsection we sketch the quasiclassical formalism used for the evaluation of the MGS properties. 
Consider an interface with normal, n, pointing in the positive x direction (n • x > 0). The incident angle, 9, of an 
electronic trajectory is defined through the relation 

kp cos 9 = • n. 

The d-wave order parameter is aligned with the crystal a-b axes according to A (A^ — kf) where k a = k^p-a, kf, = kpb. 
Introducing the misorientation angle a : n a = cos a and n b = sin a we can write k a = cos(a — 0) and k b — sin(a — 9). 
The order parameter can then be written as a function of the two angles a, 9: 

A(0) = A cos2(a-<9). (28) 

Assuming specular reflection, the momentum parallel to the interface is conserved upon reflection, while the per- 
pendicular momentum is inverted, (or equivalently, the reflected angle is given by ir — 6). The quasiclassical wave 
functions have the form of linear combinations of plane waves 

0(r,k FN ) = -^e*" 1 "" FWe***-**, (29) 



a=± 



where r = in + r 1 1 , k^ = trki?^ + k F ||, and S is the interface area. The slowly varying envelopes, (jf{x), satisfy the 
quasiclassical BdG equations, 

[w^ F -u{-id x )a z + /^e^a + + /^e-^aJ\r{x) = Er{x), j = X< ° , (30) 

where the shorthand notation is introduced, AJ = A cos 2 (<j9 — aj) , and a.j=L,R denote the angles of the a-b crystal 
axes to the normal of the interface, K = 1. It is convenient to incorporate the sign of the order parameter into the 
phase, 

3 \-¥>/2-e(-A£)7r, 3=R ^ ' 

The local interface potential is replaced in the quasiclassical approximation with the boundary conditions for slow 
wave functions envelops, 

(t> a {ir) =T° a '(t> a '(0 + ), (32) 

where T a(T denotes a general single channel (i.e. for given trajectory) transfer matrix, characterized by the transmis- 
sion amplitude, d(9), and reflection amplitude r(9), (|r| 2 + \d\ 2 = 1). 

MGS spectrum and transition matrix elements 



Here we derive equations defining the MGS energy spectrum and transition matrix elements for planar junctions 
with specular interfaces. 



The bound state solutions to equation ( 30 ) have the general form 



where 



and 7 CT = urf. + /2 with 



1 fe~^i 



, C° = J sin 77? 
J v?, ■ n J 



acos(E/|Afl), j = L 
-acos(25/|A£|), j=R 



j = L, R, 



(33) 



(34) 



(35) 



Using the properties of the spinors, the matching condition Eq. (32) can be rewritten into two sets of equations 



A° = M aa B° , = M aa B a where M aa = T aa cos (j£ - 7^ J , M aa = T aa ism (j£ - t£ J, determining the 
coefficients, A a , and, B a , upto a normalization constant. The condition that these equations have non-trivial solutions, 
detvW = 0, defines the spectral equation, 



J] sin [7g(25) - 7g(i5)] = 2? J] sin [ 7 £(25) - 7j T (25)] 



(36) 



where R 



1 -D. 



To obtain an expression for the transition matrix elements, An, we make use of the relationship with the current 
matrix elements in Eq. (18), 

where the current matrix elements are given by Eq. (|20[) within the quasiclassical approximation, 



h2 = (v£-n)(&W 



<7=± 



z=0 



2 (v£ • n)[^]*23 2 CT cos [^(250 - „£(25 2 )] 



(37) 



o-=± 



Selection rule 

For a junction with \j\ orientation, there exists a symmetry relation, A£ = A^ = Ail, and consequently, 



-Tr = Vi 



(l/2)acos(25/A a). The spectral equation then simplifies, 

cos 2 (2?7ii) = 2)cos 2 (<^/2), 



(38) 



and has the two solutions E\ = — 25 2 = v2)Aa cos(y>/2). For a spatially symmetric potential, the amplitudes in 
Eq. (33) are given by equations (upto normalization factor ensuring that (4>,4>) = 1), 



B+ 2 = VRcos((p/2), 



A+ 2 = ±y/Rcos{ip/2) 



B^ 2 = sin(277 f (25 M ) + p/2), A^ 2 = ± sm(2r/ f (25 1)2 ) + <p/2) 



(39) 



Inserting these amplitudes into equation (37) for the current matrix element (which is proportional to A), we arrive 
at the important result, 

A = 0. 

Now we show that this result is a particular case of a general selection rule forbidding transitions among the MGS 



for any symmetric junction. This selection rule is imposed by the symmetry of the Hamiltonian, Eq. (15), under 
charge and parity conjugation CV, 

CV${v) = 0*(-r). 
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To prove our statement we first note that the CP-symmetry splits the Hilbert space of the Hamiltonian (15 1, into two 



subspaces which correspond to the even and odd transformations of the eigen states under CP conjugation, 

Then we find that the operator that defines the transition matrix elements, A = iAn — — (</>i, [d ip + (i/4 : )sign(x)a z ]4>2), 
respects the CP-symmetry, and therefore the matrix elements between the states belonging to the odd-subspace and 
even-subspace vanish. 

Next, we notice that an arbitrary symmetric junction is obtained by continuous rotation of the \j\ junction. 
Such a rotation preserves the CP-symmetry, and the eigen functions transform smoothly under the rotation, unless 
the nodes of the order parameter are crossed. Therefore, the wave functions initially belonging to different discrete 



subspaces of the symmetry operator will maintain this property during the rotation. Inspection of Eqs. ( 33 1 and ( 39 ) 
for the j/j junction proves that indeed the two MGS eigen functions obtain opposite signs under CV transformation, 
and thus belong to different subspaces of the symmetry operator. 

This proves that the transition matrix elements will equal zero for the MGS of all symmetric junctions. 

J + n/j — k orientations 

The antisymmetric orientation is one of the few orientations for which one can obtain non-trivial analytical solutions 
A. For these orientations the symmetry holds, 

|A±| HASINA*. (40) 

For trajectories that admit a pair of MGS we find the spectral equation, 

cos 2 (ry + + rf) = Dcos 2 if/2. (41) 

where the shorthand 77 ± = 77^ = — 77^ was introduced for notational convenience. Using the definitions, cos 2ry ± = 
E/\A ± \, and sin2r7 ± = y/l - E 2 /\A ± \ 2 , we find the solution, 



a2 1 



, . j_ . , / — <y5 2a/1 — Dcos 2 .. 

E = ± A + A~h/Dcos- v z (42) 

1 1 2 ^(|A+| + |A-|) 2 - 4IA+A-I.DCOS 2 § V ' 

Once an analytical expression for the spectrum has been found one can also obtain an analytic expression for A in 
terms of 77 , 

_ /- (sin 27?+ - sin 27?-) / | A+A~ | sin 277+ sin 2r£ \ 
12 ' 2S sin (77+ +77-) ^ I A+ J sin 277+ + I A" J sin 277- J ' 1 ' 



Here 77* = (1/2) arccos(£'/| A* |) with the energy given by Eq. (42|. 

Notice that for k — 0, we have 77+ = 77^, so that the matrix element vanishes for this orientation, as also shown 
above. For small misorientation, k«1, this expression can be expanded into 

y/\ — D cos 2 f A 

where SA(8) = A + (0) — Aosin2# m 2kAoCOs20. This equation reduces at small transparency, 

A X i ~ VDcos^ 6 ^ w 2nv r Dcot(20). (45) 

Transition temperature 

The decay of the persistent Josephson current at large bias current applied to the junction is represented by escape 
of a fictitious particle representing the junction from a metastable potential well of the "washboard potential" formed 
by the periodic Josephson potential, Uj(ip), and potential of the current bias, U ex t((p) = —(I e /2e)ip. Grabert and 
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Weiss [2] devised a method for direct calculation of a critical temperature of transition from the thermally activated 
escape to the escape via MQT by analyzing small fluctuations around the saddle point located at the barrier top, pt- 
The semiclassical euclidian action is expanded to second order in the deviation, Sip — tp — (pt,, 



o(2) _ 



dT / dT'5ip(T)A(T-T')5<p(T r ) 

Jo 

= A(iu n )\tp n \ 2 . 



(46) 



The Fourier components, A(iv n ), of the fluctuation kernel, A(r — r'), with Matsubara frequencies, v n = 2irn/f3, 
are then the eigenvalues associated with the gaussian fluctuations around the stationary point, ipi,. In the thermal 
activation regime, the stationary point is stable, and all the eigenvalues are positive, A(iio n ) > 0. Transition to the 
MQT regime is manifested by the instability, indicated by the sign change of the smallest eigenvalue, A(ivi) < 0. The 
transition temperatures can thus be obtained by finding solutions to the equation, A(ivi) = 0. 

To evaluate the fluctuation part of the action for our system we expand in Sip and keep only second order terms. 
This gives us, 



S$\s<p]= f dT 
Jo 



c_ 

8? 



5^ 



1 d 2 U - 



2 dp 2 



where 



K(t - t') 



Sp 2 




drdr '5<p(t)K{t - t')8p>{t'), 


Vb 




) 


^SplnG- 1 





2 d(p(T)6ip{T f ) 



(47) 



(48) 



It is convenient to perform the functional differentiation in a basis where the dependence on ip is removed. This is 
achieved by a using rotation matrix U : d v U — iAU, 



K(r-r') = iTr U(r)0(r) j S(r - r') + Tr ( d(r. r'l^lr'KVi ' n^-(r) 



dp 



OH , 
dp 



(49) 



The contribution from the first part combines with the second derivative of the external potential to define the barrier 
frequency, 



8e^ 



1 d 2 U e: 

2 dp> 2 



1 



d 2 H 



dp 2 



leaving the second part which we denote by small k(r — t' 

1 

(2e) 



k(r-T>) = j^Tr (G(t,t')Ij(t')G(t',t)Ij(t) 



(50) 



(51) 



Here Ij is the current operator as defined in Eq. (64), the imaginary time (Matsubara) Green function for constant 
phase (</?(t) = pb) is given by equation, 



(Go)«(r,r') = - [6(t t')(1 p*) - 9(t> - r)p°] e"^— '>£. 



(52) 



where p® — np(E®), and E® — Ei(ipb). Due to the boundary condition, Sp(0) — 5ip(/3), we can perform integration 
by parts to obtain, 



S ( E ) [5p\ = I dTdT'Stp{ T )A(T -t')S(p(t'), 



where 



(-d 2 -U 2 )S(T-Tl) + —k(T-T>) 



(53) 



(54) 
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The Fourier representation of the operator A(r — r') defines the eigenvalues, 



H il/ n) = ^ \-(iv n ) 2 - ut - ivnloiivn)} , v n = , (55) 

8e z 1 J p 

it form, 

7o(^n) = 7^ L ^ IT " ( 56 ) 



C r / >9 9 2nn 

A(W„J = 

where — i^ / n 7o(* I/ n) = (8e 2 /C)k(iv n ) has the explicit form, 

Here £jj = -B? — £7? and ^4jj = Aij(<pb)- Comparing this result with Eq. (46) we find the equation for the transition 
temperature to the MQT regime, 

A(ifi) oc v\ — lu 2 — ifi 70(^1) = 0. (57) 

MGS and reentrance effect 

The reentrance effect described in the article results from a strong temperature dependence of dissipation produced 



by the MGS, which decrease with increasing temperature. After truncating to the MGS subspace, we present Eqs. (56) 



and (57) on the form, dropping the subscript, 

v 2 - ^1 - ii/jo(ivi) = 

2 _ 4e 2 S ^ I ieA 2 p z 



± 

8e 2 S I sA 2 p z 

v 2 1 1 



= * V V ~ C £f\±E-iv/ ~ b (58) 



C 

Here S is the junction area, and p z — np{E\) — ^(£2) while the average (. . .) is defined as 

d 2 k F , s 

(*#- (59) 

where integration is performed over the Fermi wave vectors in the positive direction of the interface normal. Assuming 
the interface to be orthogonal to the crystal a-b plane, and taking into account strong anisotropy of the Fermi surface, 
we write the integral on the form, 

,7T/2 

S(...)=N dd , (60) 

J-n/2 

where N = SIcf/Skc is the number of conducting channels for a stack 2D planes with spatial period c. For the sake of 
simplicity, we consider almost symmetric MGS spectrum, E\ = —E2 = e/2, giving p z {T) = tanh(e/4T), and proceed 
to integration over e in the integral over 8, The equation for the crossover temperature then becomes, 

where g(e) — Add/de is the MGS spectral density. The important qualitative features of the integral, independent 
of junction geometry, are the saturation effect due to the MGS population number, tanh(e/4T), and the resonance 
feature in the denominator, e 2 + (2irT) 2 . Numerical studies show that these features define the temperature dependence 
of the integral, while the role of the function, g(e)eA 2 (e), which contains information about junction geometry, is 
qualitatively insignificant. This observation allows us to approximate the latter with some constant whose magnitude 
is set by A 2 ~ D, because g(e) ~ l/e m > an d £ ~ e m ; 

r n/2 

g{e)sA 2 (e) = a / dffD(0), (62) 

J -it/2 

where a is a geometry dependent numerical constant of order ~ 1. We are then able to formulate a simple model 
equation defining the transition temperature, 



R n C Jo e 2 + (2ttT) 2 
where R n = ir/e 2 S(D) is the normal junction resistance. 
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Fitting MQT transition temperatures 

Here we shall outline the method used to fit the transition temperatures in our model to the experiment in [25] . For 



this procedure, Eq. (63 1 will be our model. Before proceeding we first simplify our notation in Eq. (63) by writing 
7o(^) = ivr\f{v '/E m ), with f(x) — f Q dytanh(iry/2x)(x 2 + y 2 )^ 1 , where r\ = 8air / R n Ce m is the coupling strength. 
The next, crucial step is to consider a dimensionless function, F(x,rf) = x 2 {l + r)f(x)), where x = 2irT/e m , and 
choose the scaling parameter e m such that the three argument values, corresponding to given transition temperatures, 
Ti, T2, and T3 give the same function value, F(xi,rj) = F(x2,tj) — F(xs,rj)] this can only be achieved by adjusting 
simultaneously the shape of the function F(x, rj) by tuning parameter 77. This procedure gives unique values for both 
parameters. Then the barrier frequency is determined by equating, uj 2 — e 2 n F(xi, rj). 

It should be noted that in our analysis we have neglected the temperature dependence of the adiabatic Josephson 
potential. In general the saturation of the MGS may lead to strong temperature dependence of the Josephson current 
- a feature suggested in [9] to be the origin of the hump structure. The model was that the potential barrier height 
changes between two asymptotically temperature independent values over a narrow region 100 mK < T < 150 mK, 
assumed to still be in the thermally activated regime. The temperature dispersion of the decay rate corresponding to 
the two different barrier heights is indicated in their Fig. 2b by two shifted lines. This explanation was consistent with 
an MQT crossover temperature T* = 50 mK obtained from the plasma frequency with estimated junction capacitance 
C = 1 pF. Later experimcnts[25 , however, suggested that this value of the junction capacitance was overestimated due 
to the presence of a stray capacitance originating from the STO substrate. Comparison with typical grain boundary 
junctions would suggest a junction capacitance of the order of 100 fF, thus increasing the crossover temperature to 
values right around the anomalous features of the temperature dispersion of the decay rate. Therefore, while the 
mechanism suggested in [S] could produce a feature like the one observed in their experiments, the parameters of this 
particular junction suggest that the reentrance effect discussed in the present paper may be more relevant. 

In addition to the stray capacitance from the substrate it was argued 25J that the c-axis transport in the tilted 
junction may cause a stray inductance. We can include the effect of such a stray LC oscillator in our analysis by 
adding an extra term, Xx 2 /(x 2 + uIq), to the function F(x,rj), where uiq — huio/e m is the dimensionless frequency of 
the stray oscillator, and A = h 2 / LqCe^ is the coupling containing the stray inductance and (unknown) capacitance 
C of the junction. The latter is connected to the barrier frequency through the relations, Wf, = uj p (1 — (Ie/Ic) 2 ) 1 ^ 4 , 
and ujp = 2eIc/frC, and evaluated through an iteration procedure, assuming switching current I e ~ 0.91c = 1.26 /xA. 
Including the LC oscillator does not produce any qualitative changes but rather slightly modifies numerical values 
of the fitting parameters. The best fit is eventually achieved for the values, e m w 320 mK, cji, « 1.7 K, lu p m 2.5 
K, and C ~ 36fF, assuming experimental values of the critical current, Ic = 1.4/iA, and the switching current, 
I e f=a 0.9/(7- Given the experimental junction transparency, D ~ 10~ 4 , we are able to evaluate the maximum energy 
gap at the interface, Ao ~ 16 K. The geometrical constant in the equation for r\ is estimated for the experimental 
value R n — 500 51, a ss 1.5, as expected. 

Linear response 

In this section we investigate the different processes contributing to the linear damping in d-wave Josephson junc- 
tions. We start by presenting the expression for the linear response in a general form, and then evaluate the contri- 
bution to the linear dissipation coming from the MGS to MGS transitions and compare that with the contributions 
coming from competing processes of nodal to nodal state transitions, and MGS to nodal state transitions. 

The non-adiabatic, real-time dynamics in Josephson junctions is described by the dynamical equations governing 
the evolution of the superconducting phase, 



C 

Ye' 



y + lf{ V ) + Tr (lj(p - p )) = I e (t), Ij = 2e (d v E + i A]) , (64) 

and the single quasiparticle density matrix, 

id t p=[H,p], H = E — ipA. (65) 



In Eq. (64), we have subtracted the adiabatic component of the Josephson current, Ij d ((fi) — Tr(Ijp°), where p° 
denotes the initial density matrix, and added an external current, I e (t) of the biasing circuit. 

The effects discussed in the article concern small oscillations around a stationary point ipQ. Straightforward lin- 



earization of Eqs. (64) and (65) with respect to small deviations from the equilibrium, (p — tpo, p — p° lead to the 
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dispersion equation, 

{-u 2 +ul+urf {u))<p u = 0, (66) 
where 70 (w) denotes the linear response of the quasiparticlcs, 

C ^TT* E»j — (w + iO) 

The indices i, j, here refer to continuous and discrete sets of quantum numbers characterizing the different eigen states, 
cj)i, of the microscopic Hamiltonian. 

MGS to MGS transitions 



Here we shall make use of the general expression for the linear response, Eq. (67), to evaluate the MGS contribution 



for nearly symmetric junctions kCI. The contribution due to MGS transitions has the form, 

4e 2 N^ r< 2 ,„ sA 2 pl 
C ^ J_ 7I /2 ±£-u-l0 

where e = Ex((po) — E 2 (tpo), pi = P11 ~ P22 ~ tanh(e/2T) and A = L4.i2(<£o)- At resonance, u < e m , where 
e m ~ \/i5Ao is the MGS bandwidth, we find the dissipative part of the linear response 

7o (u) = Imjo = (^^j ug{u)A 2 tanh ~, (69) 

where A = A((fo,9), 9 being the resonant angle, e(9) = w. To get a rough estimate of the dissipation, we use 
estimates, A ~ kVD, and <?(w) ~ l/e m , to obtain, 

///x k 2 f U! \ , ,,. 

*>~^(d ta,lh 4T- (70) 

It is instructive to compare this result to the dissipative effects of the competing processes, namely the nodal to 
nodal state transitions, and the MGS to nodal state transitions. Below we will show that both the processes produce 
dissipation, estimated at zero temperature with 

Inod ~ Inod-MGS ~ q ( ~^^J ' 

This dissipation is generally small compared to the MGS caused dissipation, by the factor ~ e m /A , unless the 
junction geometry is particularly close to the symmetric orientation. 

Nodal to nodal state transitions 



Applying the general expression, Eq ( 67 1 , for the linear response to the nodal quasiparticles we note that the 
scattering states are four-fold degenerate, and the indices i label, besides the quasiclassical trajectory angle 9 and the 
energy of the incoming quasiparticle \E\ > |A(0)|, also scattering states for electron/hole- like quasiparticles impinging 
onto the interface from the left/right labeled with s = 1, 2, 3, 4. Remembering that the temporal variation of the phase 
conserves the trajectory angle 9, we write contribution of the nodal quasiparticles to the response at zero temperature 
on the form, 

'***to--c-J_ w/a dB J J 1^— ^"(wn £ -(c+ i0 ) ' (72) 

where e = E — E' , and v(E) is the density of states, 



(73) 
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The matrix clement for the nodal state transitions is conveniently evaluated through the relation to the current, 
similar to the calculation for the MGS, 

AAE^^^^l, (74) 



I ss , (E, E') =eJ2(v F - n)C f {x, E)<j>° a , (x, E') 



(75) 



The sum over scattering states in Eqs. (72l-(75l consists of the products of various scattered waves. Focusing on 
the most interesting tunnel limit, D -C 1, we find that the main contribution is coming from the combination of the 
transmitted and reflected waves. The current for this combination, / ~ e|v> • n\y r DR, is large compared, e.g., to 
the contribution of the transmitted waves, I ~ e|vp ■ n| D. Therefore A 2 e ~ |v_p • n\ 2 RD/e, and we may present the 
matrix element factor in Eq. ( 72 ) on the form, 



^ .■> , ,s DWf ■ n| 2 , ( EE 1 
J2A 2 ss ,(E,E')e= 1 F £ 1 h [<Po,-£,-£ 



0(D 2 



(76) 



where f± is a dimcnsionlcss function constructed with the quasiclassical BdG amplitudes. With this expression the 
dissipative part of the response becomes, 



7nodM 



4e 2 N r /2 



C 



d9\A\ 



-tt/2 



dE v(E)v(Cj -E) — /i(y> , E,Cu- E). 



(77) 



where we introduced notations, E = E / A, and oj = ui/A. We see from this equation, that the resonant transitions 
at oj <C Ao select a small energy interval, 1 < E < ui — 1. This imposes a constraint on the angles, 2|A(0)| < u, 
which are restricted to small areas around the nodes. This is the source of small value of the dissipation by nodal 
quasiparticles. 

Linearizing the order parameter around this point, A « 2A #, and changing the variable, we finally obtain, 



7nodM = 



4e 2 iV oj 



duj 



dE v(E)v(u)-E)Dfi(ip ,E,L>-E). 



(78) 



C A 

This integral converges, and gives the estimate for the magnitude of dissipation produced by the nodal quasiparticles, 

1 OJ 



7nodM 



R n C A 



(79) 



MGS to nodal state transitions 

The contribution of these processes to the response function at zero temperature is given by equation, 

. 4e 2 N r' 2 f ^ g. -Aj MGS (E, E MGS )e [signjE) - sign(E MGS )} 

7nod-MGs(w) = / d6 / dEv{E) (80) 

C J e-oj-iO 

where e — E — -EmgSi and Emgs ~ VDA is the energy of the midgap state. To the lowest order in transparency, 
the components of the MGS wave functions are proportional to the factor, \J~Aj\vp ■ n|, originating from the wave 
functions normalization, and they do not depend on transparency. Evaluating the overlap in the expression for 
the current matrix elements at the transmitted side of the scattering state we find that it is proportional to e\vp ■ 
n\yD \fAj\vp ■ n| to lowest order in transparency. The transition matrix element can then be written on the form, 
similar to the nodal to nodal transitions, 

Y,<MGs(E,E UGS )e = D^I_A h (V , |, ^p) + 0(D 2 ), (81) 
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where fa is a dimensionless function of order one. Using the density of states in Eq. (73 1 and the resonance condition 
e = E — -Emgs = w we present the dissipation on the form, 

ie 2 N r' 2 jn „ fT , , IAIN {E MGS +co\ „ |A| / £mgs+w £mgs\ , co x 

J2 <Po, x >-^r~ • i 82 i 



7„od-MGs( w ) -Q- J ^d6e(E MQS +u-\A\)v^ JD 



Ul 



A 



We notice that again the resonance condition selects a small angle interval around the nodes. Furthermore, the MGS 
energy is small, -Bmgs/A ~ v~D -C 1, and can be dropped from the arguments of v, &, and fa. Then we get equation 
qualitatively analogous to the dissipation of nodal quasiparticles, 

or 



7nod-MGsM = —FT- IT I 7a "ft) h(Po,u,0), (83) 



7nod-MGsM ~ JTT^IT- ( 84 ) 



Nonlinear MGS response 

In this section we give a derivation of the nonlinear resonant response of the superconducting phase and MGS to 
the harmonic temporal oscillation of the current bias, I e (t) = I e cos ut. The starting point is the dynamical Eqs. 



(64) and (65). In what follows, we only focus on the major nonlinear effect of the MGS transitions, neglecting 
transitions between the nodal states. The MGS dynamics is described with a continuum set of two-level density 
matrices, parameterized with trajectory angle 9, and satisfying the dynamical Bloch equation, 



f>+ = (-ie - T 2 )p+ + 2tpAp z 

p z = -^A{p+ + p-) - r 1 ( Pz - pi). 



(85) 



Here the notations are introduced, p z = pn — p22, p+ = (P-)* = P12, £ = E\(ip) — £2 (</?), and the transition matrix 
element is written on the form A12 = iA(ip). Phenomenological decay rates r^Fa are introduced to account for 
intrinsic relaxation and dephasing of the MGS, due to e.g. weak short range disorder. 

We consider small oscillations of the phase around the equilibrium value ifQ, driven by the external current, I e (t) — 
I e coswi, at a frequency not far from the resonant frequency S — lj — u> p <C 1. To separate the slow and fast dynamics 
we parametrize the phase as: 

V {t) = \{^{t)e-^ + c.c.) 

2 (86) 
¥>(*) = ~Mt)e~ iut - c.c), 

where the complex variable <p u (t) = r(t)e 1 ^^ depends on the amplitude of oscillations, r(t), and the time dependent 
phase shift, &(t). Using a similar separation for the slow- and fast parts of the off-diagonal elements of the density 
matrix, 

p + (t)=p u (t)e- iut , (87) 

we get, after expanding to first order in (p — tpo and averaging over fast variables (note Ao — A(ipo) and £0 = £{<Po))] 

Pu, = -i(eo - w - iT 2 )pu, - iLL>A ip iLl p z 

■ w 1 * * \ -r ( 0\ ( 88 ) 

Pz = ^{fuPu " <P u Pu>) - T l{Pz - P z )- 

We consider the regime with slow variation of the phase oscillation envelopes, (p^, on the the time scale of the MGS 
decoherence, dt<p <C Tiip. Then the density matrix will adiabatically follow (in the rotating frame) the evolution of 
the phase amplitude, and we consider the quasi-stationary solutions, p w ,p z ~ 0, 



(uMorfGVro 



The nonequilibrium correction to the Josephson current has the form, 



Tr 



(Up 



Po, 



= 2eS(d v e{p z - p° z ) + Ae{p+ + p~) 
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(90) 

(91) 



In the linear approximation with respect to the phase amplitude r, only the last term in Eq. (91) plays a role, while 
the corrections to the diagonal matrix elements is of the second order, p z — p ~ 0(r 2 )). Expansion of this term 
recovers the result of the linear response calculation Eq. (|68|), 



(2e/C)Tr (lj(p - p )j w uj^ (uj)ip u 

with 7o(w) now containing a finite resonance broadening, r 2 , 

4e 2 S* / ujA 2 e p° z 



e lult +c.c, 



W7o(w) = 



(92) 



(93) 



C \( eo _ w )-ir 2l 

To go beyond the linear approximation we define, in analogy with the linear response analysis, a non-linear response 



(94) 



coefficient, (2e/C)Tr Ij(p - p ) ~ u-y(u, r)^^ + c.c, 



W7(w,r) 



4e 2 5 /^j ^§£o^(r) 

-c-\V°(^)-P,)+ (eo _ h;) _ ira 
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Performing integration assuming the resonance to be narrow, T2,Aot^r <C e m , we get for 7 = 7' + 47": 

7 w, r ) fa 7o = , 

"1 ,/(r^ow) 2 + r 2 

VMo^) 2 + r2' 

This is to be inserted into the averaged equation for slow phase amplitude, 

2e I 

- i2uj p Lp LJ + [-2(J P 6 + uj p j(uj p ,r)} ^ = 
-28u) p , or in a more convenient form 



(95) 
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(97) 



where I e = (e/Cw p )I e . Equations for r(t), i?(t) are obtained by dividing this equation by ip u and noticing the relation, 
tPu/^u = r/r + v&. Identifying real and imaginary parts yields the set of equations, 



r 7 "(j") ; smd 
r 2 2r 

2 2r 

The fix points of this set of nonlinear equations is found by solving equation 



-S + -7(r) 



Taking the absolute square of this equation and solving for 5 one finds, 
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(100) 



The two solutions correspond to the stable/unstable branches of the response curve. The maximum amplitude 
corresponds to the degenerate point, I 2 = (f) 2 ^, solution of this equation reads, 



T 2 a>y-p e Al 
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(101) 



